Thermal and mechanical properties of a DNA model with solvation barrier 
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We study the thermal and mechanical behavior of DNA denaturation in the frame of the mesoscopic Peyrard- 
Bishop-Dauxois model with the inclusion of solvent interaction. By analysing the melting transition of a homo- 
geneous A-T sequence, we are able to set suitable values of the parameters of the model and study the formation 
and stability of bubbles in the system. Then, we focus on the case of the P5 promoter sequence and use the 
Principal Component Analysis of the trajectories to extract the main information on the dynamical behaviour of 
the system. We find that this analysis method gives an excellent agreement with previous biological results. 

PACS numbers: 87.15.H-,87.15.A-,05.10.-a 



I. INTRODUCTION 

In the last years there exist an increasing interest in the de- 
scription of complex biological systems using simple phys- 
ical models jlj. Physical models for biological problems 
should contain the key ingredients to explain their basic phe- 
nomenology . In this sense, DNA melting can be modeled 
using very simple statistical mechanics models J2]. One of 
the most successful is the one dimensional Peyrard-Bishop- 
Dauxois (PBD) model UHH which with very few assumptions 
is able to reproduce the melting curves for different DNA se- 
quences. PBD model undergoes an entropic phase transition 
between the native closed state and a denatured one in which 
the DNA strands are separated. The nature of this transition 
strongly depends on model parameter, so different versions of 
the model have been proposed in the last years. 

A virtue of the PBD model is that it allows to study not only 
the equilibrium properties of the molecule but also dynamical 
ones; for instance, the formation and stability of the so called 
DNA bubbles (short open segments of the DNA chain) JHH]. 
Due to its versatility, this simple model has been applied in 
other contexts beyond the study of the equilibrium melting 
curves. Thus, it has been used to investigate the DNA me- 
chanical denaturation J3,Ht] and it has been specially success- 
ful in the modelling of the open regions in short DNA hair- 
pins J(| |H [T^l ■ This fact has motivated to extend the model 
to study non-homogeneous sequences and to try to correlate 
the dynamical behavior of the model to functional aspects of 
the DNA chain. To be specific, it has been suggested that 
the localised dynamical excitations (bubbles) of the model are 
directly related with the protein-DNA binding sites flU [l2tl . 
This point has been controversial, due to the simplicity of the 
PBD model assumptions EH. 

One of the most interesting improvements to the model 
is the inclusion of a barrier in the intra-base interaction 
term 1 1 311 - Such barrier accounts for solvent interactions in 
the system. The addition of barriers modifies deeply the na- 
ture of the denaturation transition and the dynamics of the 
molecule. Depending on the parameter values, the transition 
become sharper, even in the harmonic stacking energy ver- 
sion. The transition width has not a trivial behavior and it 



deserves a more careful study. On the dynamical side, the 
bubbles lifetime increases dramatically with the inclusion of 
this term and approaches to the experimental observations. 

The purpose of this paper is twofold. First, we character- 
ize the melting transition of homogeneous sequences and the 
formation of bubbles phenomenon for a wide range of pa- 
rameters (section IV). This allows us to set suitable param- 
eter values for next research. Then, we focus on the problem 
of the relation between dynamics and function. In section V, 
we study a heterogeneous sequence which contains known in- 
formation on the transcription process of DNA to mRNA (a 
P5 virus promoter). We use Principal Component Analysis 
of trajectories to identify the regions in the sequence which 
most relevantly contribute to the dynamical fluctuations of the 
molecule. These regions are the softest ones, from a mechan- 
ical point of view, and they correlate fairly well with the rel- 
evant biological sites. Given its simplicity and efficiency, we 
propose to extend this analysis to other sequences in order to 
get insight over functionality of the different genome regions. 



II. THE MODEL 

We study a modification of the Peyrard-Bishop-Dauxois 
model. In the PBD model the complexity of DNA is re- 
duced to the study of the dynamics of the N base pairs of the 
molecule. For each base pair we define the variable y n associ- 
ated to the distance between the bases. The total energy of the 
system is then approached by: 
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2m 



+ V(y n )+W(y n ,y n - 1 ) 



(1) 



Here p n = m dy„/dt, n is the index of a base pair and m its 
reduced mass. 

In this equation, we identify two energy potential terms: 
V(y„), an on-site potential one, and W(y„,y„-i) which ac- 
counts for inter-pairs interactions. The standard PBD model 
corresponds to a particular choice for these energy terms. 

The potential V(y„) describes the interaction between the 
two bases of a pair. The PBD model uses a Morse potential to 
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FIG. 1: (Color online) Morse potential with (solid lines) and with- 
out (dashed lines) solvation barrier for A-T (black) and C-G (or- 
ange/light gray) base pairs. For A-T pair without barrier a = 4 A 
and D = 0.043eV. For A-T with barrier a = 4A _1 , D = 0.05185eV, 
G = 3D, yo = 2/ a and b = 1 /2a 2 (the value of D has been adjusted 
in each case to give the desired melting temperature for the A-T uni- 
form chain). For C-G pairs Dqq = 1.5£>at and occg = 1-5ccat- 



account for such interaction (see Fig. Q]). This is an standard 
approximation for chemical bonds. 



The melting process involves the breaking of hydrogen 
bonds between the bases. In order to study different DNA se- 
quences the PBD model also can include sequence dependent 
Morse potential parameters: D„ and oc„. It can be intuitively 
deduced from the fact that A-T pairs are linked by two hy- 
drogen bonds meanwhile C-G ones are linked by three, thus 
forming a more stable link. Following fl6ll in our simulations 
we will use Dcg = 1-5Dat an d OCcg = 1.5(Xat- Although re- 
cent works have considered more complex sequence depen- 
dence we will not introduce more complexity at this 
level in the model. 

We will work with a modified version of the PBD model to 
include solvent interaction. Solvent interaction stabilizes open 
pair states by means of the hydrogen bonds that the base pairs 
may form with the solvent when opened lfl5ll . Such bonds 
have to be broken before the pair closes again. This effect can 
be included in a simple way with an effective barrier in the 
Morse potential. The addition of the barrier has to avoid, as 
far as possible, any other effect in the Morse well. Hence, we 
have chosen the following definition for the on-site potential 
V(y): 



V(y)=D(e- ay ~l) 2 + Ge- 



-yoV/b 



(4) 



V(y)=D(e 



l 



(2) 



This potential results from the addition of a gaussian barrier, 
whose height is controlled by G, the position is given by yo 
and its width by b. A reasonable election for such parameters 
is G = 3D, yo = 2/oc and b = 1 /2a 2 . Figure [T]plots the intra- 
base potential V for an A-T and a C-G base pair. 



where D corresponds to the dissociation energy of the pair 
[V(0) = and V(°°) = D] and a sets the amplitude of the po- 
tential well [V"(Q) = 2D(X ], See also that variable y measures 
deviations of the base distance with respect to equilibrium, 
thus V(y) has a minimum for y = 0. 

An important contribution to the stability of DNA molecule 
comes from the stacking of the bases. The potential 
W(y„,y„-i) describes the interaction between base pairs along 
the DNA strand. The simplest approximation is to con- 
sider just harmonic interaction between adjacent pairs in the 
molecule. However this choice does not reproduce the obser- 
vation of a sharp thermal denaturation curve. To solve this, the 
PBD model includes a position dependent nonlinear coupling 
constant 

W(y n ,y n -i) = ^(l + pe- 5 ^-')) (y.^y,^) 2 . (3) 

The effect of this term, whose intensity is governed by p is 
to change the effective coupling constant from K(l + p) to K 
when one of the base pairs is displaced away from its equilib- 
rium position. The parameter 8 sets the scale length for this 
behavior. 

DNA melting or denaturation refers to the separation of 
the two strands of the DNA molecule to generate two single 
strands. The great merit of the PDB model defined in Eqs. (O 
and (O is to reproduce the different melting curves observed 
in DNA molecules. Thus, it has been used as a solid starting 
point for many other statistical mechanics studies of DNA. 



III. METHODS 

A. Langevin Dynamic Simulations. 

In order to study the behavior of the system we have 
performed molecular dynamics numerical simulations of the 
Langevin equation: 
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d\W(yn,yn-i)+W(y n+h y„)] 
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Ut): (5) 



where m is the mass of the pair, y the effective damping of the 
system and cj(?) accounts for thermal noise, (^„(f)) = and 
(^„(f)^t(f')) = 2myk B Tb nk b(t -t') with T the bath tempera- 
ture. 

The equations were numerically i nteg rated using the 
stochastic Runge-Kutta algorithm um , [1911 . Simulation of 
large A-T chain (section IV) used periodic boundary condi- 
tions in order to avoid any terminal effect. For the P5 pro- 
moter (section V) we used fixed boundary conditions which 
will be described below. 

To best characterize the phase transition we have computed 
the mean energy (u), and mean displacement (y) defined as 



j N,t s 

< M > = 77T £[WW»-i) + VCy,,)] 
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Nt, 



(7) 



is the number of base pairs in the sequence to study and t s 
the total simulation time. 

We are also interested in computing the probability of being 
opened for each base. We define this as follows: 



fl.(yni) = 7L®(y»(0-yai) 
fj t 



(8) 



where y t h is a threshold for opening that we choose just behind 
the barrier of the Morse potential and &(x) is the Heaviside 
step function (&(x) = for x < and &(x) = 1 for x > 0). 
Thus P n indicates the fraction of time the n-th base is opened 
during the time t s . 



B. Principal Component Analysis (PCA) 

PCA is a statistical method to extract information from a 
large set of data in a multidimensional phase space, allowing 
to reduce the dimensionality of the variables to those that in- 
clude most of the fluctuations of the original system |20ll . This 
is achieved by a change of coordinates in the phase space and 
the new coordinates are the so-called Principal Components, 
PCs. This method has become a standard tool in the analysis 
of trajectories in "all-atom" simulations of macromolecules 
& 

Operationally, we have to build the N x N correlation ma- 
trix 



C{i,j) = (yiyj)-(yi)(yj)- 



(9) 



Diagonalising this matrix we obtain an ordered set of eigen- 
values (A-i > A-2 > A.3 . . .) with their corresponding eigenvectors 
(vi, V2, V3....). Eigenvalue A.,- gives the amount of fluctuations 
which corresponds to eigenvector v,, thus the new coordinates 
are ordered in such a way that the few first ones retains most 
of the fluctuations of the system. It is useful to define a fre- 
quency CO, associated to each eigenvalue X, and given by 



CO, 




(10) 



From this point of view, the larger fluctuations will correspond 
to the lowest frequencies. It can be proved that in the low tem- 
perature limit these frequencies converge to the normal mode 
frequencies of the system. Normal mode analysis has also 
used in the detection of relevant motions in coarse-grained 
model of proteins l22tl . 



IV. THE UNIFORM CHAIN 

A. Fitting the phase transition: parameters of the model 

In this section we will study the melting transition of a ho- 
mogeneous chain of A-T pair bases. Although there are some 
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FIG. 2: (Color online) The melting transition. Mean energy and dis- 
placement (inset) as a function of temperature of a homogeneous 
chain for the standard (open circles) PBD model and the modified 
one (full circles). Parameters are given in the text. 



studies on the influence of the different parameters in the tran- 
sition, up to our knowledge, there is not a systematic scan 
over the parameter space. In particular, when the barrier is in- 
cluded on the Morse potential a clear narrowing of the transi- 
tion is observed (see Fig. [2]). The same effect is expected when 
the p parameter is increased. Thus, one of the first problems 
when dealing with this modified PBD model is to find suitable 
model parameters. This is mainly done by adjusting the the- 
oretical melting temperatures to the experimental ones. The 
melting temperature T,„ is usually defined as the temperature 
at which the DNA strands are half denatured. For a homo- 
geneous A-T chain, the melting transition has been reported 
around T m = 310A" Q10. 

One should be aware that, in general, the melting tem- 
perature depends on the length of the sequence, base com- 
position, topological structure and salt concentration. Thus, 
it is difficult to predict the exact transition temperature of a 
given sequence. In the case of the uniform chain the melting 
temperature depends strongly on salt concentration, which is 
not explicitly included in this mesoscopic model. Usual pa- 
rameters for the standard PBD model can be found in ref- 
erences iflrjl \\% for instance. Following previous findings 
and guided by our numerical simulation results we choose 
the following set of values for studying the standard model: 
to = 300 Da, D = 0.043eV, a = 4A~\ K = 0.03eVxA 2 , 
p = 3 and 8 = 0.8A (these values are slightly different 
to those given in lfl7lo . If barrier is included we modify 
the value of D to obtain the same T m . Thus, in this case 
D = 0.05185eV, and G~3D~ 0.1556eV, y = 2/oc = 0.5A 

and b = 1 /2a 2 =0.03 125 A 2 . 

The melting transition curve for the model with and without 
solvation barrier is shown in Fig. [2] for a homogeneous A-T 
sequence with N = 220. At high temperature is observed the 
expected linear behavior of a free gaussian polymer chain; i.e. 
the completely unzipped state. 
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FIG. 3: (Color online) Transition temperature versus the stacking 
constant K for various p values. 



FIG. 5: (Color online) PC frequency spectrum at different temper- 
atures between 120 and 330K. Frequency units are (Z)/m) 1 / 2 a = 
5.15 x 10 12 sec -1 . 
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FIG. 4: (Color online) Transition width versus the stacking constant 
K for various p values. 
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The numerical determination of the melting temperature is 
not a trivial issue. For the case of a uniform chain, we use 
the following computational criteria to determine the melting 
temperature and transition width from the (u(T)) and (y(T)) 
curves. First we define two temperatures. The larger one, T2, 
gives the onset of the linear behavior in (u(T)) which indi- 
cates that the chain is completely melted. The other one, T\ 
estimates the beginning of the transition defined as the tem- 
perature for which (y(T)) = yo, i.e. when the chain can be 
considered to be by average in the barrier position. Then we 
define the transition width AT = Tj_ — T\ and the melting tem- 
perature as T m = (T\ +T2) /2. 

It is interesting to briefly discuss on the effect of the main 
model parameters in the melting transition and the dynamical 
behavior of the system. High values of D, the Morse poten- 
tial dissociation energy, produces also high T m . Parameters 
K and p, set the stacking interaction between adjacent bases. 
Its effect in T,„ and AT can be seen in Figs. [3] and [4] A suit- 



FIG. 6: (Color online) Temperature dependence of the lowest PC 
frequency. The straight line corresponds to the best fit with the 
parameters given in the text. Frequency units are (D/m)^ 2 a = 
5.15 x 10 12 sec _1 . 



able melting temperature and transition width is obtained at 
high K values (thus we chose K = 0.03eV/A 2 ) and moderate p 
(p ~ 3). Moderate p are biologically satisfactory and avoid the 
very narrow bubbles observed for high p values. With respect 
to the solvation barrier its presence makes bubbles live longer. 
Thus a complete separation of the strands is also facilitated 
leading to a decrease of the melting temperature (this change 
can be counterbalanced by an increasing in the D value). The 
inclusion of the barrier also reduces the width of the transition, 
allowing the election of moderate p values. 
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FIG. 7: Typical molecular dynamics simulation trajectories for the PBD model of a homogeneous sequence without barrier (upper figures) 
and with barrier (lower ones). Dark areas correspond to open base pairs. Long-living bubbles are clearly observed when solvation barrier is 
included. Small figures show a base pair configuration at a given time (left) and the trajectory of a given base as a function of time (right). 
Trajectory time is 200 ns and y„ is given in units of a -1 = 0.25A. 



B. PCA of the phase transition 

PCA provide us another method to analyse the phase transi- 
tion. Figure shows the spectrum of effective frequencies as 
defined in Eq. ( fTOb for the uniform A-T chain. At low temper- 
atures, a sharp band of linear waves excitations close to the 
ground state is clearly identified. At very low tem pera tures, 
the frequencies are given by the dispersion relation 112411 : 

mor ~2Da 2 + 2K(l +p)[l -cos(nn/N)]. (11) 

Thus, the gap for the lowest mode is controlled by the intra- 
base potential parameters. In this limit the PC eigenvectors are 



the normal modes of a homogeneous array. As temperature 
increases, the nonlinear excitations are more important. Such 
nonlinear modes are responsible of the larger fluctuations and 
correspond to the larger eigenvalues or lower frequencies. As 
we approach the melting temperature a soft mode goes to zero 
(see Fig. |5), the chain is not pinned by an on-site potential 
and every strand moves freely. Hence at high temperature the 
frequencies are those of a free gaussian chain with coupling 
given by the stacking parameter K: 



mco 2 =2K[1 -cos(nn/N)]. (12) 



6 



Note that the zero mode corresponds to the largest wave- 
length. 

Figure [6] shows the temperature evolution of this mode fre- 
quency. The curve can be fitted using a critical behavior func- 
tion CO (T m - T) v with T m = 0.51 and v = 2.2. The char- 
acterization of the order of the phase transition is a difficult 
issue l25l l26tl . To our knowledge dynamic exponents for this 
family of models are not known. 



C. Bubble formation 

One of the more interesting aspects of DNA molecule 
amenable to be studied in the framework of the PBD mode, 
is the formation and stability of the so called DNA bubbles 
(short open segments of the DNA chain). The onset of such 
states have an important role in the understanding of the oper- 
ation of DNA molecules. 

We observe that the inclusion of a solvation barrier in the 
model has a dramatic effect in the dynamics of the system 
increasing the lifetime of bubbles. Figure|7]compares molec- 
ular dynamics numerical simulations for a model with barrier 
to another without barrier and at p = 3. As expected the pres- 
ence of barriers modifies the kinetics of base opening and fur- 
ther closing. The individual trajectories of the different bases 
clearly show this behavior. Although the profile of the bub- 
bles is similar in the two cases, the non-barrier dynamics is 
completely different from that with barrier. Without a barrier 
the opening of the base pair corresponds to a large amplitude 
oscillation along the Morse potential, and an easy closing is 
favoured. With barrier, the kinetics is controlled by the pres- 
ence of two equilibrium states separated by the solvation bar- 
rier. Closing events are more difficult and bubbles live longer. 
Although in this way, we approach to experimental values of 
lifetime bubbles (of the order of ns), this fact makes the sim- 
ulations more difficult since it is required very long runs to 
have good statistics. 

We see that parameter p affects the cooperativity between 
open bases and then the bubbles lifetime and width. A large 
p produces long live bubbles but extremely narrow ones (one 
or two bases), as those shown in reference ll27tl with p = 25. 
With this value of p, the transition is extremely narrow and, 
after a long transient, bubbles nucleate to drive the unzipping 
of the whole chain. A moderate value, p = 3, is good enough 
to get both longer and wider bubbles (four to ten bases), as 
shown in figures. 



V. THE P5 PROMOTER SEQUENCE 

A promoter is a region of the DNA sequence where the tran- 
scription to mRNA is initiated and controlled. In prokaryotic 
cells, regions of -10 and -35 bp (upstream from the transcrip- 
tion start site) are the most important sites for transcription 
regulation. In this section we analyse the P5 core prom oter 
which has been widely studied in the literature lfl2[ l28ll . We 
show that the PCA of Langevin trajectories of the molecule 
clearly identifies the biologically relevant sites. 
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FIG. 8: (Color online) PC frequency spectrum at different tem- 
peratures for the P5 promoter sequence. Frequency units are 
(D/m) 1 /2 a = 5.l5x 10 12 sec _1 . 
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FIG. 9: (Color online) Top figure: Probability (not normalised) of 
opening for sequence of P5 promoter. Bottom figure: The three first 
PC eigenvectors at T = 290K. 



The sequence of the P5 promoter is given by the 69 pb: 5'- 
GTGCCCATTTAGGGTATATATGGCCGAGTGAGCGAGC 
AGGATCTCC ATTTTGACCGCAAATTTGAACG-3 ' . We 
will perform numerical simulations of this sequence (pa- 
rameter values for C-G and A-T base pairs have been given 
previously). In order to avoid unphysical denaturation events 
at low temperature due to finite size effects |@], we apply 
the following boundary conditions: First we add 10 C-G 
bp sequences to the ends of the promoter to create hard 



time 



FIG. 10: Trajectory of the P5 sequence promoter at T = 29QK. The larger bubbles appear mainly in the region of the localised PC eigenvectors 
showed in Fig. [9] The time of the showed trajectory is 400 ns. 



boundaries. Second, we set the extremes to (closed) to 
avoid the complete opening of the chain. 

Figure[8]shows the evolution of the PC frequency spectrum 
with temperature. Several features distinguish these curves 
from the homogeneous case. For low temperatures, two bands 
corresponding to A-T (lower band) and C-G (upper band) 
links are clearly identified. The eigenvectors are localised 
around rich regions in both kinds of complementary pairs. At 
intermediate temperatures the gap between both bands disap- 
pears and the eigenvectors broaden. C-G pairs surrounded by 
A-T ones are more likely to open and frequencies diminish. 

Close to transition (T ps 345^T), several modes detach to 
low frequencies. The first three modes are strongly localised 
and show peaks in four regions (see Fig.|9]l. These modes rep- 
resent the "softest" regions of the sequence, are related with 
high probability of opening and drive the unzipping. 

This picture is validated by the study of the three first eigen- 
vectors (the large eigenvalue or small frequency ones). Fig- 
ure [9] shows the measured probability of opening in the se- 
quence, defined by Eq. ©, and the computed eigenvectors. 
An excellent correlation between both figures is observed. 
Region +1 corresponds to the transition starting site (TSS) 
whereas regions -30 and -40 correspond to the binding sites 
of transcription factors i.e. a region rich in A/T pairs like the 
TATA-box. Localised eigenvectors span over regions of ten 
base pairs which fairly correspond with the width of the bub- 
bles. 

Figure [10] plots a typical trajectory, showing a few bubbles 
mainly at the regions pointed by the PC eigenvectors. Note 
that at this temperature very few opening events occur, so a 
precise statistics on bubble formation becomes difficult. How- 
ever, PC analysis gives a good account of these sites even with 
the presence of a few bubbles. Interestingly, the third eigen- 
value shows a clear correlation between region close to TATA 
box and the TSS. This represents motion fluctuations that in- 
volve both regions. Thus, modification of one of them could 
interfere in the fluctuations in the other opening a channel for 
the regulation of transcription. 

VI. CONCLUDING REMARKS 

We have studied the dynamics of a DNA molecule using the 
PBD model with solvation interaction. Solvation interaction 



is modeled by including a barrier in the usual on-site Morse 
potential. First we have analyzed the melting transition in an 
A-T homogeneous chain. We conclude that the inclusion of 
the barrier not only modifies the phase transition but also has a 
great influence in the dynamics of localised excitations (bub- 
bles). In combination with the non linear parameter p, the 
width and lifetime of bubbles can be tuned. We find a set of 
parameters which will be suitable for future studies with this 
model. 

We have applied the model to an inhomogeneous sequence 
with biological meaning, a virus promoter. The study of tra- 
jectories using Principal Component Analysis allows us to 
detect the biological relevant sites without the need of long 
molecular dynamics runs. We have found that the softest 
modes of the PC spectrum are highly localised in those sites. 
Even more, these sites have been also detected as the most 
likely to be opened in bubbles 112911 . In this work we cannot 
elucidate the controversy on whether the DNA can direct its 
own transcription as was suggested in lfl2ll and then make use 
of theoretical methods to obtain the functional regions. How- 
ever, we can obtain from the mesoscopic PBD model some 
regions which are more sensible to the formation of bubbles 
(as extracted from the eigenvectors of the PCA). Of course, 
the simplicity of the model can not take into account other 
effects, like the flexibility of DNA molecule, but it is able to 
tackle those related with the sequence. 

Finally, we stress the use of a tool, the PCA, to the anal- 
ysis of statistical mechanics simulations of simple models to 
obtain useful information on linear and nonlinear excitations 
and the relation with its phase transition. 
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